% X \in S^{n(p+1)}
% C \in S^{n(p+1)}
% TODO decompose into functions? not sure how to do that without losing cvx scope
% TODO vectorize; this works but it slow
function B = regular_mlm(C, n, p, gamma)
	cvx_begin
		variable X(n*(p+1), n*(p+1)) symmetric
    % construct D
		expression D(n, n*(p+1));
		for i=0:p
			D(1:n,1:n) = D(1:n,1:n) + X((1:n)+n*i,(1:n)+n*i);
        end		
        for k=1:p
			for i=0:(p-k)
				D(1:n,(1:n)+n*k) = D(1:n,(1:n)+n*k) + 2*X((1:n)+n*i,(1:n)+n*(i+k));
			end
        end
    % evaluate penalty h
		expression h(1);
		for i=1:n
			for j=(i+1):n
				h = h + norm([D(i,j+n*(0:p)) D(j,i+n*(0:p))], Inf);
			end
		end

		minimize(-log_det(X(1:n,1:n)) + trace(C*X) + gamma*h);
		X == semidefinite(n*(p+1));
	cvx_end

    B0 = sqrtm(X(1:n,1:n));
    B1p = B0 \ X(1:n,(1+n):(n*(p+1)));
    B = [B0 B1p];
end
